home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / Tools / Dev / bignum / bignum2.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-29  |  6.6 KB  |  219 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /* BIGNUM2.C
  4. **
  5. ** Routines to do Big Number Arithmetic. These are multi-precision, unsigned
  6. ** natural numbers (0, 1, 2, ...). For the storage method, see the BigNum
  7. ** typedef in file BIGNUM.H
  8. **
  9. ** Released to the public domain by Clifford Rhodes, June 15, 1995, with
  10. ** no guarantees of any sort as to accuracy or fitness of use.
  11. */
  12.  
  13. #include <stdlib.h>
  14. #if defined(_MSC_VER)
  15.  #include <memory.h>
  16. #elif defined(__TURBOC__)
  17.  #include <mem.h>
  18. #else
  19.  #include <string.h>
  20. #endif
  21.  
  22. #include "bignum.h"    /* Typedef for BigNum type */
  23.  
  24. BigNum * BigNumAlloc(int nlen)
  25. {
  26.       /* Allocates memory for a BigNum object with nlen entries. Returns a
  27.       ** pointer to the memory with the data array initialized to zero if
  28.       ** successful. If not successful, NULL is returned.
  29.       */
  30.  
  31.       BigNum * big;
  32.  
  33.       big = (BigNum *) malloc(sizeof(BigNum));
  34.       if (big != NULL)
  35.       {
  36.             big->nlen = nlen;
  37.             if (nlen < 1)
  38.                   big->n = NULL;
  39.             else if ((big->n =
  40.                   (UShort *) calloc(nlen, sizeof(UShort))) == NULL)
  41.             {
  42.                   free(big);
  43.                   return NULL;
  44.             }
  45.       }
  46.       else  return NULL;
  47.  
  48.       return big;
  49. }
  50.  
  51. void BigNumFree(BigNum * b)
  52. {
  53.       /* Frees memory used by BigNum object b. */
  54.  
  55.       if (b != NULL)
  56.       {
  57.             if (b->n != NULL)
  58.                   free(b->n);
  59.             free(b);
  60.       }
  61. }
  62.  
  63. BigNum * BigNumDiv(const BigNum * a, const BigNum * b, BigNum ** c)
  64. {
  65.       /* Big Number a is divided by Big Number b. If successful a pointer to
  66.       ** the quotient is returned. The user must supply a pointer to a Big
  67.       ** Number pointer, c, to receive the remainder.
  68.       ** If unsuccessful, NULL is returned.
  69.       ** Neither a nor b is changed by the call.
  70.       */
  71.  
  72.       int      i, j, d, bpos;
  73.       UShort   carry, quo;
  74.       long     m1, m2, m3;
  75.       BigNum * quotient, * atemp, * btemp;
  76.  
  77.       /* Find non-zero MSB of b, make sure b is not 0 */
  78.  
  79.       for (bpos = 0; bpos < b->nlen && b->n[bpos] == 0; bpos++)
  80.             ;
  81.       if (bpos == b->nlen)  /* Attempt to divide by zero! */
  82.             return NULL;
  83.  
  84.       /* Create a copy of b, making sure btemp->n[0] != 0 */
  85.  
  86.       if ((btemp = BigNumAlloc(b->nlen - bpos)) == NULL)
  87.             return NULL;
  88.       memcpy(btemp->n, b->n + bpos, btemp->nlen * sizeof(UShort));
  89.  
  90.       /* Create a copy of a, making sure atemp->n[0] == 0 */
  91.  
  92.       carry = (a->n[0] == 0) ? 0 : 1;
  93.       if ((atemp = BigNumAlloc(a->nlen + carry)) == NULL)
  94.       {
  95.             BigNumFree(btemp);
  96.             return NULL;
  97.       }
  98.       memcpy(atemp->n + carry, a->n, a->nlen * sizeof(UShort));
  99.  
  100.       /* Allocate memory for quotient and remainder */
  101.  
  102.       if ((quotient = BigNumAlloc(max(1, atemp->nlen - btemp->nlen))) == NULL)
  103.       {
  104.             BigNumFree(atemp);
  105.             BigNumFree(btemp);
  106.             return NULL;
  107.       }
  108.       if ((*c = BigNumAlloc(btemp->nlen)) == NULL)
  109.       {
  110.             BigNumFree(atemp);
  111.             BigNumFree(btemp);
  112.             BigNumFree(quotient);
  113.             return NULL;
  114.       }
  115.       d = MODULUS / (btemp->n[0] + 1);
  116.       for (carry = 0, j = atemp->nlen - 1; j >= 0; j--)   
  117.       {
  118.             m1 = ((long) d * (long) *(atemp->n + j)) + (long) carry;
  119.             *(atemp->n + j) = (UShort) (m1 % (long) MODULUS);
  120.             carry = (UShort) (m1 / (long) MODULUS);
  121.       }
  122.       for (carry = 0, j = btemp->nlen - 1; j >= 0; j--)
  123.       {
  124.             m1 = ((long) d * (long) *(btemp->n + j)) + (long) carry;
  125.             *(btemp->n + j) = (UShort) (m1 % (long) MODULUS);
  126.             carry = (UShort) (m1 / (long) MODULUS);
  127.       }
  128.       for (j = 0; j < (atemp->nlen - btemp->nlen); j++)   
  129.       {
  130.             if (*btemp->n == *(atemp->n + j))
  131.                   quo = MODULUS - 1;
  132.             else
  133.             {
  134.                   m1 = (long) *(atemp->n + j) * (long) MODULUS;
  135.                   m1 = (m1 + (long) *(atemp->n + j + 1)) / (long) *btemp->n;
  136.                   quo = (UShort) m1;
  137.             }
  138.             m3 = (long) *(atemp->n + j) * (long) MODULUS +
  139.                   (long) *(atemp->n + j + 1);
  140.             do
  141.             {
  142.                   if (btemp->nlen > 1)
  143.                         m1 = (long) quo * (long) *(btemp->n + 1);
  144.                   else  m1 = 0l;
  145.                   m2 = (long) quo * (long) *btemp->n;
  146.                   m2 = (m3 - m2) * (long) MODULUS +
  147.                         (long) *(atemp->n + j + 2);
  148.                   if (m1 > m2)
  149.                         quo--;
  150.             } while (m1 > m2);
  151.  
  152.             bpos = btemp->nlen - 1;
  153.             i = j + btemp->nlen;
  154.             m2 = 0l;
  155.             while (i >= j)
  156.             {
  157.                   if (bpos >= 0)
  158.                         m1 = (long) quo * (long) *(btemp->n + bpos);
  159.                   else  m1 = 0l;
  160.                   m3 = (long) *(atemp->n + i) - m1 + m2;
  161.                   if (m3 < 0l)
  162.                   {
  163.                         m2 = m3 / (long) MODULUS;
  164.                         m3 %= (long) MODULUS;
  165.                         if (m3 < 0l)
  166.                         {
  167.                               m3 += (long) MODULUS;
  168.                               m2--;
  169.                         }
  170.                   }
  171.                   else  m2 = 0l;
  172.                   *(atemp->n + i) = (UShort) (m3);
  173.                   bpos--;
  174.                   i--;
  175.             }
  176.             if (m2 == 0l)
  177.                   *(quotient->n + j) = quo;
  178.             else
  179.             {
  180.                   *(quotient->n + j) = quo - 1;
  181.                   bpos = btemp->nlen - 1;
  182.                   i = j + btemp->nlen;
  183.                   carry = 0;
  184.                   while (i >= j)
  185.                   {
  186.                         long tmp;
  187.  
  188.                         if (bpos >= 0)
  189.                               carry += *(btemp->n + bpos);
  190.                         tmp = carry + (long) *(atemp->n + i);
  191.                         if (tmp >= (long) MODULUS)
  192.                         {
  193.                               tmp -= (long) MODULUS;
  194.                               carry = 1;
  195.                         }
  196.                         else  carry = 0;
  197.                         *(atemp->n + i) = (UShort) tmp;
  198.                         bpos--;
  199.                         i--;
  200.                   }
  201.             }
  202.       }
  203.       j = atemp->nlen - btemp->nlen;    
  204.       bpos = 0;
  205.       carry = 0;
  206.       while (j < atemp->nlen)            
  207.       {
  208.             m3 = (long) carry * (long) MODULUS + (long) *(atemp->n + j);
  209.             *((*c)->n + bpos) = (UShort) (m3 / d);
  210.             carry = (UShort) (m3 % d);
  211.             j++;
  212.             bpos++;
  213.       }
  214.       BigNumFree(atemp);   /* Free temporary memory */
  215.       BigNumFree(btemp);
  216.  
  217.       return quotient;
  218. }
  219.